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Abstract 



We present simulation results on the equilibrium relaxation of Brownian 
planar rotors based on a uniformly frustrated XY model on a square lat- 
tice. The rotational relaxation exhibits typical dynamic features of fragile 
supercooled liquids including the two-step relaxation. We observe a dynamic 
cross-over from high temperature regime with Arrhenius behavior to low tem- 
peraure regime with temerature-dependent activation energy. A consistent 
picture for the observed slow dynamics can be given in terms of caging effect 

and thermal activation across potential barriers in the energy landscapes. 
PACS No.: 64.70Pf, 05.45.+j, 64.60.Cn 
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I. INTRODUCTION 



Last decade or so have witnessed significant advances in our understanding of the under- 
lying mechanism for the slow dynamics of supercooled liquids approaching the glass transi- 
tion . The development of mode-coupling theory of supercooled liquids and extensive 
experiments and computer simulations have played crucial roles in such advances. Some 
efforts have also been devoted to devise model systems (even though somewhat artificial) 
M which show glassy behavior similar to that of supercooled liquids. One line of research 
along this direction is to find (lattice) model systems with no quenched disorder but some 
intrinsic frustration built into the model, which may exhibit glassy relaxations [|5|-§|. 

One can imagine that there may exist a common microscopic mechanism which underlies 
the observed similarities in the relaxations of model systems and real supercooled liquids. 
This possibility is made more plausible by the universal scaling property observed in the 
dielectric susceptibilities of a variety of supercooled liquids and some plastic (glassy) crys- 
tal [p!0|-[T2]|. Here in this work, we address the question of this possible common mechanism 
by investigating the equilibrium orientational relaxation of planar Brownian rotors whose 
interaction is prescribed by that of uniformly frustrated XY (UFXY) models with dense 



frustration, which is a prime example of non-randomly frustrated systems [|n| characterized 
by complex degeneracy of ground states and many metastable states. 

While a recent simulation || of the present authors deals with the relaxation of the vortex 
charge density for a purely dissipative dynamics, here we examine directly the orientational 
relaxation with finite rotational inertia, which offers more transparent views on the origin 
of the observed slow relaxation. Also, due to the one-dimensional nature of the phase of 
the planar rotors, it is convenient to probe the properties of the angular motions of the 
rotors of the system. We find that, by including phenomeno logical rotational inertia in the 
dynamic equation for the rotors, the orientational correlation exhibits a two-step relaxation, 
which is analogous to the (fast) (3 and a relaxations of supercooled liquids. Mean square 
angular displacement (MSAD) exhibits three stage behavior, i.e., the early time ballistic, 
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intermediate sub-diffusive, and late time diffusive regimes, which is argued to be consistent 
with the picture of the cage effect and long-time activated dynamics for the motion of the 
rotors. It is shown that there exist two dynamically distinct regimes: a high temperature 
regime where the dynamics is governed by a temperature-independent activation energy, 
and a low temperature regime, in which the activation energy increases with decreasing 
temperature, which is interpreted as arising from complex energy landscapes [14 , 15] probed 
by the system in the low temperature regime. 

II. DYNAMIC MODEL AND SIMULATION METHOD 

We consider the following Langevin dynamics for a collection of planar rotors on a square 
lattice 

Iu t (t) + ^(t) = + Vl (t) (1) 

where / is the moment of inertia, uJi(t) = 6i(t) the angular velocity of the rotor at site i, 7 
the damping constant, and rji (t) the thermal noise. The equation (|I|) describes the Brownian 
motion of rotors subject to the interaction potential energy ^({0}). The thermal noise 77, (t) 
is given by a gaussian random variable 

< vM > = 

< r}i{t)m{t') > = Z-yTSijSit - t') (2) 

where the Boltzmann constant ks is set equal to unity. The variance of the noise in (^|) 
ensures that the system at temperature T evolves toward the equilibrium state whose prop- 
erties are governed by the Boltzmann distribution exp(— E({8}, {u>})/T) where the energy 

E({9}, M) is given by E({9}, M) = I + V{{6}). 

Here we chose the potential energy as the energy of the two dimensional UFXY 

model on a square lattice which takes the form pll 



V({0}) = -JE cos (^ - 0; - A-,-) (3) 



where J is the coupling constant and (ij) denotes nearest neighbor pairs. The bond angles 
satisfy the constraint 

E A a = ( 4 ) 

where the sum is over belonging to the unit plaquette P causing competing interactions 
(frustration) between the rotors. Here, / is called the frustration parameter of the system. 

A convenient choice for A^ is the Landau gauge which is given by Ay = for every 
horizontal bond and Ay = ±27r/xj for the vertical bond directed upward (downward) with 
Xi being the x-coordinate of the site %. It can be readily checked that this choice of the bond 
angles obeys the condition fl4j). Due to the invarince of the Hamiltonian ([I]) under / — > / + 1 
and / — > —f, we need to consider the values of / only over the range [0, 1/2]. A physical 
realization of this model can be found in the two dimensional square array of Josephson 
junctions under a uniform perpendicular magnetic field. In this situation, the bond angle 
A^ is identified with the line integral of the vector potential A of the transverse magnetic 
field: Aij = (27r/$ ) f? A ■ dl where $ is the flux quantum <3> = hc/2e per unit plaquette. 
With this identification the strength of magnetic field B is given by Ba 2 = /$o where a is 
the lattice constant. 

The UFXY model can be mapped jL7| onto that of a lattice Coulomb gas with charges 



of magnitude (n — f), n — 0, ±1, ±2, • ■ -, where charges correspond to phase-vortices with 
suitably defined vorticity around the plaquettes. The lowest excitation consists of charges 
with magnitudes 1 — / and — /, respectively. The charge neutrality condition then implies 
that the number density of positive charges is equal to /. For the case of / = 0, the well- 



known Kosterlitz-Thouless transition |18| occurs via vortex-antivortex unbinding at a finite 
temperature. Except for this case of unfrustrated XY model, the equilibrium nature and 
associated phase transitions of these systems are not very well understood even for the next 



simplest case of / = 1/2, the so-called full frustrated XY model |19[. For example, the 
ground state configurations for the case of general / = p/q (p and q are relative primes) are 
not known [^,0] except for some low order rational values of /, such as / = 1/2, 1/3, 2/5, 



3/8, etc, where staircase type of ground state configurations are known analytically p2| , |2T 



As q becomes large (the limit of irrational frustration), due to the complexity of the 
degeneracy of the system and long equilibration time, it is quite a difficult task to analyze 
the nature of the low temperature phase of the system. And, inspite of recent claim by 
Denniston and Tang P31 that there exist a first order transition near T c ~ 0.13J, in the case 
of / = 1 — g, (g being the golden-mean ratio g = (y/E— 1)/2 ~ 0.618), it is fair to say that the 
low temperature phase is not completely understood yet. On the other hand, since it is clear 
that many metastable states are possible due to the dense frustration, one can expect that 
Brownian dynamics ([I]) with the potential energy (0) may generate a slow relaxation where 
trapping of the configurations in deep metastable minima and thermal activation across the 
potential barriers play a crucial role. Note that there is no intrinsic disorder in the present 
system, which distinguishes itself from a spin glass system where both intrinsic disorder and 



frustration are considered to be essential 24 



With the potential energy the Langevin equation is explicitly given by 

Iui(t) + iu)i(t) = -j£sin(0j - 6j - A i3 ) + m(t) (5) 
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We integrate the equation fl5|) in time, starting from random initial conditions {^(0)} and 
{co>j(0)} using an Euler algorithm on a square lattice of linear size N = 34, In our simulations, 
we used I = 1.5, 7 = 1, J = 1 and / = 13/34, which is a Fibonacci approximant to 
/ = 1—g. Periodic boundary conditions are employed for both spatial directions. The 
results were averaged over 150 ~ 1000 different random initial configurations, depending 
on the quenching temperature. As for the integration time step, we used dt = 0.05 in 
the dimensionless unit of time. No essential difference could be found in the results when 
compared with those obtained by using dt = 0.01. 



III. RESULTS AND DISCUSSIONS 



In order to probe the orientational relaxation of the system we first computed the on-site 
auto-correlation function for the planar spins 
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C R (t) = ± (^cos&iO) - Hi))) (6) 



where the bracket < • • • > in (^) represents an average over different random initial configu- 
rations. In this work we focus only on the lowest order correlation even though one may also 
measure the higher order correlations, as was done in recent molecular dynamics simulations 
2|-g7|. 



Shown in Fig. 1 is the on-site auto-correlation function Cn(t). The relaxation continu- 
ously slows down as the temperature is lowered. In order to characterize the slowing down 
of the relaxation, one can define a characteristic relaxation time Tr(T) as Cr(tr) = 1/e. 
The temperature dependence of Tr(T) is shown in the inset of Fig. 1. It exhibits an Ar- 
rhenius behavior at high temperatures, while at low temperatures (T < 0.20) it shows 
a non-Arrhenius behavior, which can be well fitted by the Vogel-Tamman-Fulcher form 
r R (T) = r exp[DT /(T - T )] with r ~ 9.92, T ~ 0.08, and D ~ 3.58 (28|. Similar 
non-Arrhenius behavior was observed in the vorticity relaxation as well [§]. 

An interesting feature of the rotational relaxation is that it exhibits a two-step relaxation, 
a very fast relaxation (up to t ~ 3 for T = 0.13 J, the lowest temperature probed) and a slow 
relaxation following the fast relaxation. The earliest part of the fast relaxation is expected 
to be well described by the free rotation of the rotors Iui(t) + 70^ (t) = 0. For the time 
range where t <C J, the inertial term is dominant and hence 6i(t) — 6i(0) ~ u^(0)£. It is then 
easy to show that the relaxation is given by Cn(t) ~ 1 — (T/2I)t 2 using the equipartition 
theorem (uj 2 ) = T/I. 

The long-time part of the slow relaxation can be well fitted by the stretched exponential 
form Cji(t) = C exp[— Ci^/tr)^] (C\ = 1 + lnC due to the definition of r R ), shown in 
Fig. 2. We find that the exponent f3 varies with temperature: it decreases as the temperature 
is lowered, as shown below in the inset of Fig. 3. It is interesting to note that at low 
temperatures (T < 0.2) the short time part of the slow relaxation shows a deviation from its 
stretched exponential fit and the time region for this deviation tends to extend over longer 
time regions with lowering temperature. We have fitted this region with a power law decay 
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known as the von-Schweider relaxation Cr(£) = C2 — Cst b where the exponent b also 
varies with temperature (see the inset of Fig. 3). We now examine the scaling behavior of the 
rotational relaxation. Shown in Fig. 3 is Cn(t) versus the rescaled time £/tr(T). Obviously 
the earliest part of the relaxation does not obey the scaling since faster time scale (the inverse 
of the inertia which is temperature independent) is involved in this regime. We also observe 
that the time-temperature superposition of the relaxation function is systematically violated 
in the late (slow) part of the relaxation, especially at low temperatures. This breakdown of 
the scaling is consistent with the fact that the two exponents b and (3 vary with temperature. 

It would be interesting to examine the response function corresponding to the orien- 
tational correlation function Cn(t). The response function in the frequency (u) domain 
can be defined as (via fluctutation dissipation theorem) x ( u ) = 2vrz/ / °° dt cos(2tt ist)Cji(t). 
Fig. 4 shows x'iy) versus v in a semi-log plot. We see that there exist two peaks, the 
low-frequency a peak and the high-frequency peak (microscopic peak). As the temperature 
is lowered, the a-peak moves to lower frequency, indicating the slowing-down of the reori- 
entational relaxation. At the same time, the maximum value of x"(^), which is analogous 
to the Debye- Waller factor, continuously decreases, and the a-spectrum becomes broadened 
as the temperature is lowered. We also note that as the temperature is lowered a minimum 
of the spectrum is slowly developed. All these features in the frequency spectrum of the 
orientational relaxation is qualitatively quite similar to the recent broad-band dielectric sus- 
ceptibility measurement of supercooled liquids P, pO| , pT| . According to the recent dielectric 
susceptibility data, the a-spectrum of supercooled liquids consists of two power law regimes 
in the right-hand side of the a-peak. The first power law relaxation clearly corresponds 
to the stretched exponential relaxation in time domain. In addition to this, another power 
law regime is observed in the high frequency side of the a-spectrum. It is quite interesting 
that similar power law relaxation is also observed in the high frequency part of the mag- 



netic susceptibility of a spin glass system [32]. Although we can not better resolve the high 



frequency part of the a-spectrum of the present orientational relaxation due to the bad 
statistics of the spectrum at low temperatures, we believe that our orientational relaxation 



spectrum also exhibits similar two-power-law regimes in the right hand side of the a- peak. 
The reason is that, even though the long time part of Cn(t) can be well fitted by a stretched 
exponential function, the regime of its validity (for stretched exponential form) is limited 
to late time regime only and does not extend to intermediate time regime where so called 
von-Schweidler relaxation (with different exponent b) better fits the relaxation function. 
In the frequency domain this will correspond to two power law behavior. 

In order to investigate the self-diffusion of the rotors, we measured the mean squared 
angular displacement (MSAD) 

((A0(i)) 2 > = ^ (jb(0i(t) - 0i{O))^ (7) 

where the phase angle 9i(t) is unbounded. Fig. 5 shows a log-log plot for the MSAD 
((A9(t)) 2 ) versus time t. For all temperature range probed, we see that ((A9(t)) 2 ) ~ t 2 
in the early time regime, which may be called the ballistic regime. It is expected that 
each rotor makes a free rotation in this time regime. Hence the MSAD is then given by 
((A9(t)) 2 ) ~ (T/I)t 2 in the ballistic regime. This regime corresponds to the earliest part 
of the relaxation Cn(t) ~ 1 — (T/2I)t 2 . For high temperatures this ballistic regime directly 
crosses over to the diffusive regime where ((A9(t)) 2 ) ~ t. But as the temperature is lowered, 
in the intermediate time regime a sub-diffusive regime characterized by ((A9(t)) 2 ) ~ v' 
with cf) < 1 (for example, <fi ~ 0.3 for T = 0.13 J) starts to appear and extends over more 
than two decades of time at the lowest temperature probed (T = 0.13 J). The sub-diffusive 
regime sets in at the same time t ~ 2 for all temperatures. In this regime the rotational 
motion is significantly hindered. This can be directly seen in Fig. 6 which shows the an- 
gular displacements A9i(t) = 9i(t) — 6^(0) at some representative sites at T = 0.15J. We 
clearly see from this figure that for all these phase angles the rotational motion looks almost 
frozen for more than a few thousand time units. This strongly indicates that the system is 
stuck in a particular configuration among many possible metastable states. The rotor then 
executes a local vibrational motion only, which corresponds to the caging in the dynamics 
of real supercooled liquids. At longer time scales, however, the local rotors can execute full 
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rotations via activated tunneling through the potential barriers, showing occasional abrupt 
rotational motions, as shown in Fig. 6. Similar jump motions have been observed in MD 
simulations of soft-sphere mixtures | 53fl , binary Lennard- Jones p3| , and the colloidal glass 



Also, neighboring rotors can execute collective rotations, thereby slowly rearranging 
the whole phase configurations. This stage will correspond to the slow part of C R (t). This 
entire time evolution of the self rotational motion is qualtitatively the same as that observed 
in MD simulations of the orientational relaxation of molecular supercooled liquids |27 . 



The rotational diffusion constant Dr(T) can be obtained by the slope of the MS AD ver- 
sus t in the long time limit where MSAD exhibits diffusive behavior ((A6(t)) 2 ) = 2D R (T)t. 
As shown in Fig. 7, at high temperatures the rotational diffusion constant exhibits an Ar- 
rhenius behavior, which is well fitted by D R (T) = D exp(—AE/T) with D ~ 0.68 and the 
temperature independent activation energy AE ~ 0.87J. As the temperature is lowered, 
however, Dr(T) shows a strong deviation from the Arrhenius behavior. This behavior im- 
plies that the long time dynamics in the high temperature regime is governed by activation 
barriers whose average height does not depend on temperature. In the low temperature 
regime, the rotors explore deeper valleys in the potential energy landscapes whose depth 
increases as the temperature decreases, giving rise to the non-Arrhenius behavior of the 
relaxation time |37 . 



It was observed in some experiments of supercooled liquids |38] that while both transla 



tional and rotational diffusion constants are proportional to the inverse of viscosity at high 
temperatures, the decrease of the translational diffusion constant is less dramatic than the 
inverse of viscosity at low temperatures. The rotational diffusion constant, on the other 
hand, is still proportional to the inverse of viscosity at low temperatures down to the glass 
transition. This relative enhancement of the translational self-diffusion is also revealed in re- 
cent simulations of supercooled liquids [§9|,^0j and the lattice model systems pTT| , |42|| . Here we 



compared the temperature dependences of the two time scales 1/Dr(T) and tr(T). Shown 
in the inset of Fig. 7 is a plot for D r {T)tr{T) versus T. Since the product D r {T)tr(T) 
in the plot is measured to be nearly contant down to T = 0.20J, the two time scales are 
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observed to be proportional to each other, i.e., tr(T) ~ D^iT)^ 1 up to T = 0.20 J. The 
data points below 0.20J tend to deviate from this proportionality, indicating more rapid 
decrease (rather than enhancement) of the rotational diffusion constant. However, it is not 
clear to us whether this anomalous behavior is a genuine feature of the present model or 
not. 

We have also measured the normalized angular velocity auto-correlation function (AVCF) 



CMit) = EEW(o» ■ (8) 



In the absence of the interaction between rotors, Cav{1) can be easily obtained as Cav(^) = 
exp(— 'jt/T). With interaction, as shown in Fig. 8, the AVCF shows a strongly damped 
oscillatory motion. As the temperature is lowered, the amplitude of oscillation becomes 
enhanced. This behavior strongly indicates that the rotors execute angular rattlings in 
'cages' |3[. 



For purely gaussian distribution of the angular displacements, it is easy to show that the 
rotational correlation function Cn(t) can be expressed in terms of the mean square angular 
displacement ((A8(t)) 2 ) as cj^(t) = exp(-((A#(£)) 2 )/2). Shown in Fig. 9 is the comparison 
of the rotational correlation function Cn(t) and its gaussian approximation Cjf\t). We find 
that Cr{€) exhibits a good agreement with the gaussian approximation in the early time 
regime whereas it shows a considerable deviation from the gaussian approximation in the 
late time regime. In order to characterize the non-gaussian nature of the distribution of 
displacements, the non-gaussian parameter has often been used in simulations of supercooled 
liquids [fPj -fr7|. Here we measure the same quantity for the angular displacements, which is 
defined as 

= i {mm _ , (9) 

V ; 3((A^(t)) 2 ) 2 V ; 

where the factor 1/3 comes from the one dimensional nature for the motion of the rotors. 
As shown in Fig. 10, a^i) exhibits three time regimes of distinct behavior, as in the MSAD. 
It almost vanishes in the ballistic regime and then rapidly increases toward its maximum 
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in the intermediate time regime, and finally decreases again in the long time regime. This 



temporal behavior is qualitatively the same as that observed in some MD simulations E6 



As the temperature is lowered, the maximum value of ot2{t) rapidly increases, and at 
the same time, the time regime where atiit) increases are extended, indicating strong non- 
gaussian nature of the rotational motion in this regime. This regime corresponds to the 
sub-diffusive regime in the time dependence of the MS AD shown in Fig. 5. It is expected 
that a2(t) eventually decays to zero since, for pure diffusion, the gaussian distribution is 
expected for the angular displacement. 



IV. SUMMARY 

We have shown that the relaxation of a phenomenological Brownian rotors based on 
densely frustrated XY model Hamiltonian exhibits a slow dynamics which is remarkably sim- 
ilar to the relaxation of fragile supercooled liquids. We find that there exist a dynamic cross- 
over from high temperature regime where the dynamics can be described by temperature- 
independent activation energy, and low temperature regime where non-Arrhenius behavior 
sets in, which can be attributed to the dynamic characteristics of the system probing deeper 
valleys in the potential energy landscapes with increasing height of the activation energy 
barrier. The caging in the metastable minima and thermal activation across potential barri- 
ers in the energy landscapes may provide the underlying physical origin for the similarity in 
the slow dynamic behavior of the present model system and that of real fragile supercooled 
liquids. It would be very interesting to quantitatively characterize the metastable states 
present in the system such as finding the local minima and densities of metastable states. In 
this regard, it would also be very instructive to examine how the dynamic features change as 
the value of the frustration parameter / is varied. We can also consider Newtonian dynam- 
ics version of our system and compare with Langevin dynamics p£8||i9|| , which may provide 
further insight into these questions. We will undertake further study along these directions 
in the near future. 
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FIGURE CAPTIONS 



Fig. 1. The rotational auto-correlation functions C R (t) versus time t (in dimensionless units 
with 7 = 1 and J = 1 ) for temperatures T/J = 0.5, 0.4, 0.3, 0.25, 0.2, 0.17, 0.15, 
0.14, 0.13. Inset: An Arrhenius plot for the characteristic relaxation time defined as 
C(t r (T)) = 1/e, where solid line is a Vogel-Tamman-Fulcher fit at low temperature 
regime (see the text). 

Fig. 2. Stretched exponential fit (dashed lines) to the long time part of the autocorrelation 
functions (for the same temperatures as in Fig. 1). Time t is measured in the same 
dimensionless units as in Fig. 1. 

Fig. 3. Rotational autocorrelation functions C R (t) versus the rescaled time t/r R (T). Note 
that the time-temperature superposition is systematically violated. The inset shows 
the temperature dependence of the exponents b(T) and (3(T) characterizing the slow 
part of the correlation function C R (t). 

Fig. 4. Dynamic response function x" ( u ) corresponding to the rotational relaxation versus 
frequency v for temperatures T = 0.5, 0.4, 0.3, 0.25, 0.2, 0.17, 0.15. In addition 
to the microscopic peak, one can clearly see the development of /5-minimum (as the 
temperature is lowered), decrease of the height of the a peak and broadening of the 
width of the a peak. 

Fig. 5. Mean squared angular displacement ((A6(t)) 2 ) versus time t (in dimensionless units) 
for the same temperatures as in Fig. 1. At the lowest temperature probed (T = 0.13 J), 
sub-diffusive regime extends over more than two decades. 

Fig. 6. Angular displacement A0j(i) versus time t (in dimensionless units) at some chosen 
lattice sites for T = 0.15 J. Rotational caging effect and occasional jump motions are 
exhibited. 



18 



Fig. 7. An Arrhenius plot for the rotational diffusion constant D R (T). We can see a crossover 
from high temperature regime with Arrhenius behavior to low temperature regime with 
non-Arrhenius behavior. The inset shows an anomalous deviation from the Stokes- 
Einstein relation by plotting the product D r (T)tr(T) versus T, where we can find 
that, at low temperaures, the coefficient of angular diffusion is smaller than that which 
would be expected from standard Stokes-Einstein relation. 

Fig. 8. The angular velocity auto-correlation functions CAvif) for T = 0.50 J and T = 0.13 J (t 
in dimensionless units). For comparison, dotted line represents exponential relaxation 
corresponding to the situation where the potentials are neglected. One can see a strong 
rotational cage effect indicated by the oscillating tail of CAvif)- 

Fig. 9. The rotational autocorrelation functions versus time t (in dimensionless units) for tem- 
peratures T/J = 0.5, 0.3, 0.17, 0.14, and 0.13 together with Gaussian approximation 
results (dotted lines). Systematic deviations are seen at late time stage. 

Fig. 10. Nongaussian parameter versus time t (in dimensionless units) for the same tempera- 
tures as in Fig. 1. 
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